function corrp = osl_clustertf(c,thresh,nP,con)

% corrp = osl_clustertf(c,thresh,nP,TFCE)
%
% c is data where:
% nS = size(c,1); %number of 'subjects' (recording sites)
% nF = size(c,2); %number of frequencies
% nT = size(c,3); %number of timebins
%
% thresh is threshold to use on c
%
% nP is num of permutations
%
% con is the connectivity criterion. Could be 6(surface)
%               18(edge) or 26(corner). For the 2D case here, these will
%               correspond to 4, 8 and 8 respectively.
%
% LH and MWW 2012

c=permute(c,[1 4 2 3]);

if nargin<2
  thresh = 3.1;
end

if nargin<3
  nP = 500;
end

TFCE = 0; %TFCE support - not yet implemented

nS = size(c,1); %number of 'subjects' (recording sites)
nR = size(c,2); %number of regressors
nF = size(c,3); %number of frequencies
nT = size(c,4); %number of timebins

cr = reshape(c,nS,nR*nF*nT);

%build up a null distribution of the maximum cluster size for each
%permuatation i and each regressor r:
for i = 1:nP
  dm = ((rand(nS,1)>0.5)-0.5)*2;
  
  [cg,vg,tg] = ols(cr,dm,eye(size(dm,2)));

  tg = reshape(tg,size(dm,2),nR,nF,nT);

  for r = 1:size(tg,2)
    [imlabel LL] = spm_bwlabel(double(squeeze(tg(1,r,:,:))>thresh),26);
    tmp = unique(imlabel);
    tmp(tmp==0) = [];
    nL = 0;
    if ~isempty(tmp)
      for k = tmp' %loop over clusters
	nL(k) = sum(sum(imlabel==k));
      end
    end
    nulldist(i,r)=max(nL);
  end
end

%run a one sample t-test on data and compare cluster size to
%permutation data
dm = ones(nS,1);
[cg,vg,tg] = ols(cr,dm,eye(size(dm,2)));

tg = reshape(tg,size(dm,2),nR,nF,nT);

corrp = zeros(size(tg));
for r = 1:size(tg,2)
  [imlabel LL] = spm_bwlabel(double(squeeze(tg(1,r,:,:))>thresh),26);
  tmp = unique(imlabel);
  tmp(tmp==0) = [];
  cpimlabel = zeros(size(imlabel));
  if ~isempty(tmp)
    for k= tmp' %loop over clusters
      nL = sum(sum(imlabel==k));
      cp = mean(nL>nulldist(:,r));
      cpimlabel(imlabel==k)=cp;
    end
  end
  corrp(1,r,:,:) = cpimlabel;
end